#ifndef SimpleMultiDimensionGaussian_HPP
#define SimpleMultiDimensionGaussian_HPP
#include "matrixdecomposition.hpp"
#include "randomNumberGenerators.hpp"
#include <vector>
#include <cmath>

namespace TermStructure
{
    //! A Gaussian Generator based on CholeskyDecomposition
	template<typename TVarianceMatrix,typename T,typename TSize = long>
	class SimpleMultiDimensionGaussian
	{
        //! A Gaussian number generator
		GaussGenerator<T> _gaussian;
        //! the storage for the generated value
		mutable std::vector<T> GaussianVector;
        //! Covariance matrix decomposed by Cholesky's method
		SimpleMatrix<T,TSize>  _cholDecompositionMatrix;
	public:

		//! Set Variance Matrix
		void SetVarianceMatrix(const TVarianceMatrix &_varianceMatrixinput)
		{			
			_cholDecompositionMatrix.Resize(_varianceMatrixinput.GetLength(0),_varianceMatrixinput.GetLength(1));
			CholeskyDecomposeTo(_varianceMatrixinput,_cholDecompositionMatrix);
			GaussianVector.resize(_varianceMatrixinput.GetLength(0),0);
		}
		//! Generate Next Random Vector
		void NextGroup() const
		{
			for(TSize i=0;i<_cholDecompositionMatrix.GetLength(0);i++)
			{
				GaussianVector[i]=_gaussian.Next();
			}
			for(TSize j=_cholDecompositionMatrix.GetLength(0)-1;j>=0;j--)
			{
				GaussianVector[j]=GaussianVector[j]*_cholDecompositionMatrix(j,j);
				for(TSize k = j-1; k>=0;k--)
				{
					GaussianVector[j]+=GaussianVector[k]*_cholDecompositionMatrix(j,k);
				}
			}

		}
		//! Get the i-th element of the generated vector
		const T &operator()(TSize i,TSize j) const
		{
			return GaussianVector[i];
		}
		//! Get the size of the generated vector
		const TSize GetLength(int i) const
		{
			return this->GaussianVector->size();
		}

	};
    //! A Reduced-rank Gaussian Generator utilizing SVD. Rank is defaul to 4
	template<typename TVarianceMatrix, typename T, typename TSize = long>
	class PCABasedMultiDimensionGaussian
	{
        //! the desired rank. If the size of the variance matrix is smaller, the size of the variance matrix is used
		TSize _rank;
        //! U of the SVD, A=UDV^T
		SimpleMatrix<T,TSize> u;
        //! squared of D of the SVD, A=UDV^T
		SimpleMatrix<T,TSize> squareRootOfd;
        //! The generated random vectors
		mutable std::vector<T> GaussianVector;
		mutable std::vector<T> _temp;
        //! a gaussian number generator
		GaussGenerator<T> _gaussian;
	public:
        //! Constructors. Rank is defalt to 4
		PCABasedMultiDimensionGaussian(TSize rank = 4)
		{
			SetRank(rank);
		}
        //! Set the rank of the 
		void SetRank(TSize rank)
		{
			this->_rank = rank;
		}
        //! Set Variance Matrix
		void SetVarianceMatrix(TVarianceMatrix &_varianceMatrixinput)
		{
			u.Resize(_varianceMatrixinput.GetLength(0),_varianceMatrixinput.GetLength(1));
			this->squareRootOfd.Resize(_varianceMatrixinput.GetLength(0),_varianceMatrixinput.GetLength(1));
			SimpleMatrix<T,TSize> d(_varianceMatrixinput.GetLength(0),_varianceMatrixinput.GetLength(1));
			SimpleMatrix<T,TSize> v(_varianceMatrixinput.GetLength(0),_varianceMatrixinput.GetLength(1));
			SingularValueDecomposeTo(_varianceMatrixinput,u,d,v);

			for (TSize i=0;i<(d.GetLength(0)<_rank?d.GetLength(0):_rank);i++)
			{
				squareRootOfd(i,i)=sqrt(d(i,i));
			}
			_temp.resize(d.GetLength(0)<_rank?d.GetLength(0):_rank);
			GaussianVector.resize(squareRootOfd.GetLength(0));
		}
		//! Generate Next Random Vector
		void NextGroup() const
		{
			for(TSize i=0;i<(this->squareRootOfd.GetLength(0)<_rank?squareRootOfd.GetLength(0):_rank);i++)
			{
				_temp[i]=_gaussian.Next();
			}
			for(TSize j=0; j<squareRootOfd.GetLength(0);j++)
			{
				GaussianVector[j] = 0;
				for(TSize i=0;i<(squareRootOfd.GetLength(0)<_rank?squareRootOfd.GetLength(0):_rank);i++)
				{
					GaussianVector[j]+=u(j,i)*squareRootOfd(i,i)*_temp[i];
				}
			}
		}
		//! Get the i-th element of the generated vector
		const T &operator()(TSize i,TSize j) const
		{
			return GaussianVector[i];
		}
        //! Get the size of the generated vector
		const TSize GetLength(int i) const
		{
			return this->GaussianVector->size();
		}
	};
}
#endif